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ABSTRACT 

The evolution of a pulsar wind nebula (PWN) inside a supernova remnant (SNR) is sensitive to 
properties of the central neutron star, pulsar wind, progenitor supernova, and interstellar medium. 
These properties are both difficult to measure directly and critical for understanding the formation of 
neutron stars and their interaction with the surrounding medium. In this paper, we determine these 
properties for PWN G54.1+0.3 by fitting its observed properties with a model for the dynamical and 
radiative evolution of a PWN inside an SNR. Our modeling suggests that the progenitor of G54.1+0.3 
was an isolated ~ 15 — 20 M 0 star which exploded inside a massive star cluster, creating a neutron star 
initially spinning with period Pq ~ 30 — 80 ms. We also find that > 99.9% of the pulsar’s rotational 
energy is injected into the PWN as relativistic electrons and positrons whose energy spectrum is well 
characterized by a broken power-law. Lastly, we propose future observations which can both test the 
validity of this model and better determine the properties of this source - in particular, its distance 
and the initial spin period of the central pulsar. 

Subject headings: pulsars: individual: PSR J1930+1852, ISM: individual objects: PWN G54.1+0.3, 
ISM: supernova remnants, X-rays: individual: PWN G54.1+0.3 


1. INTRODUCTION 

Stars born with a mass > 8 M 0 (e.g., iHeeer et al.l 
120031 ) are believed to end their lives in a core-collapse su¬ 
pernova powered by the gravitational c ollaps e of its iron 
core into a neutron star (e.g. IZwickvlll938f) . In many 
cases, this collapse creates a rapidly spinning (initial ro¬ 
tational period Pq <C 1 s) neutron star with a strong 
(B ~ 10 12 G) surface magnetic field observed as a pul¬ 
sar. The rotational energy of such neutron stars powers 
a mag netized , highly relativisti c outflow called a pulsar 
wind (IGoldreich fc Julianl 119691 : 1 Aronsl 120021) . The con- 
finements of this outflow c reate s a “termination shock” 
(IKennel fc Coronitil Il984bt see IGaensler fc Slanel 120061 
for a recent review), and the post-shock (“downstream”) 
pulsar wind creates a pulsar wind nebula (PWN) as it 
expands into its surroundings. When the neutron star 
is very young, it is located inside the supernova rem¬ 
nant (SNR) created by the expansion of the material 


the properties of the progenitor supernova (e.g., the 
mass and initial kinetic energy of the superno va ejecta; 
IKennel fc Coronitil 1 1984al I Gelfand et al.ll2009l and refer¬ 
ences therein) - quantities difficult to measure directly 
but vital for understanding the physics of core-collapse 
supernovae. 

Currently, the best way of measuring the proper¬ 
ties of the central neutron star, its pulsar wind, and 
progenitor supernova requires modeling the dynami¬ 
cal an d radiative evolution of a PWN inside an SNR 
(e.g. [Reynolds fc ChevalieU 11984 IGelfand et all 120091 : 


iTanaka fc Takaharal 20ldT Bucciantini et al.l 1201 1L again 
see IGaensler fc Slanel 120061 for a recent review!. Such 
models have been developed, incorporating the effect 
of the spin-down of the central neutron star (e.g. 
iBucciantini et al.l 12004 IGelfand et al.l 120071 1200911 , the 
evolution o f the surrounding SNR as it expands into the 
ISM (e.g. IGelfand et al.l [2007112009H . and for different 
nronerties of the nulsar winds after beinff iniected into 










































possible parameter space and id entify degeneracies be¬ 
tween parameters, and in Section 13.2.11 compare our de¬ 
rived properties of the neutron star, pulsar wind, progen¬ 
itor supernova, and surrounding ISM with the results of 
previous analyses to determine the impact of our assump¬ 
tions. In Section [4] we discuss the implications of these 
results concerning the progenitor of this system (Section 
14.111 . the formation of its associated pulsar (Section |4.2L 
and both the production and acceleration of particles in 
the pulsar wind ('Section [4.3D . Finally, in Section [5J we 
use our model to predict the results of future observa¬ 
tions of this source and discuss their potential implica¬ 
tions. Lastly, in Section [6l we summarize our results. 


2. EVOLUTIONARY MODEL 

Our model for dynamical and radiative evolution of a 
P WN inside an SNR is closely based on that developed 
bv Idelfand et al.i (2009). We assume the rotational lu¬ 
minosity E of the cent ral neutron star evolves as (e.g. 
iGaensler fc Slanell2006f) 


E(t ) — Eq 



(1) 


where t is the time since the progenitor supernova, Eq 
is the neutron star’s initial spin-down luminosity, r s d is 
neutron star’s “spin-dow n” timesca l e, an d p is the neu¬ 
tron star’s braking index dGelfand et al.l|20 09), and that 
all of the rotational energy of the neutron star is carried 
away by the pulsar wind generated in its magnetosphere. 
We further assume that, immediately after the pulsar 
wind is injected into the PWN at the termination shock, 
a constant fraction 7 /b of its energy is in the form of mag¬ 
netic fields, while the rest 1 — Pn is in the kinetic energy 
of electrons and positrons (|Gelfand et al.l 12009t h Theo¬ 
retical studies predict that, under most physical condi¬ 
tions, the spectrum of these particles is well described by 
a relativis tic Maxwellian w i th a high-energy power-law 
tail (e.g., iSpitkovskvl 120081 : ISironi fe SoitkovsTvl 1 201 ll) . 
While this spectrum can reproduce the broadband spec¬ 
tral ene rgy di s tribu tion (SED) of some PWNe (e.g., 
iFang fo Zhanel 12010H . we find it does not work for 
G54.1+0.3 for constant parameters. Instead, we use a 
broken po wer-law inject spectrum, which recent simula- 
tions (e.g.. ISironi fe Si)itkovskvl[20 1 It ISironi et alJl2013lf 
are able to produce under certain physical condition, and 
has been used to reproduce the broadb and SE D of this 
PWN and others in similar work (e.g. [Chevalier] 120051 
iBucciantini et al.ll2011b FTbrres et al .1120 l ll l . In this case, 
the injected particle spectrum is: 


calculate n b reak(t) by requiring that: 

(1 — 77 B )E= J En(E)dE (3) 

-£'min 

at all times. 

To minimize the number of free parameters, we as¬ 
sume that all parameters related to the properties of the 
pulsar wind (?7b, .Emin, Ebreak, E max , pi, p 2 ) remain con¬ 
stant with time - in contrast to other models which as¬ 
sume different temporal ev olution’s for (som e ) of t hese 
parameters. For example, IBucciantini et all ( 201 1! as¬ 
sume that E max is proportional to the electric potential 
of the pulsar’s magnetosphere <h, while others set E max to 
the particle energy whose Larmor radius is the equal to 
the th e radius of the terminati on shock (e.g-.lTorres et all 
120141) or the PWN itself (e.g.- lLi et al.ll20 1 01) . In Section 
13.2.11 we estimate the systematic uncertainty resulting 
from these different assumptions by comparing our re¬ 
sult to those derived from different models. We also note 
that our model does not consider the possibility of ions 
in the pulsar wind, nor the magnetic reconnection and 
particle acceleration beyond (“downstream”) of the ter¬ 
mination shock as pre dicted by recent 3D sim ulations 
of these systems (e.g., iPorth et all 12013112014D . These 
processes are expected to p rimarily affect th e spectral 
evolution of the PWN (e.g., lOlmi et aD[2014 1. and are 
left for future work . 

As done bv iGelfand et all (I2009H . the dynamical evo¬ 
lution of the PWN is determined by the motion of the 
surrounding shell of swept-up material. This shell is sub¬ 
ject to a net force resulting for the difference in pressure 
between the PWN and the SNR just outside the PWN. 
We calculate the pressure just outside the PWN using 
the procedure described bv IGelfand et al.l (12009H . which 
assumes the initial density profile of the supernova ejecta 
is a uniform density inner core surrounded by an enve¬ 
lope whose density decreases as p oc r~ 9 , where r is the 
distance from the center of the SNR, and that the SNR 
is expanding in to a constant dens ity ISM. 

As done bv lGelfand et al.l (12009!) , we calculate the pres¬ 
sure inside the PWN assuming that both the PWN’s 
magnetic field strength E pwn and the particle density are 
spatially uniform i.e., using a “one-zone” model for the 
PWN. We account for both adiabatic and radiative losses 
(assumed to be dominated by synchrotron emission and 
inverse Compton scattering off Cosmic Microwave Back¬ 
ground photons) of the electrons and positrons inside the 
PWN. The spectrum of photons generated by the radia¬ 
tive losse s are calc ulated using the same procedure de¬ 
scribed in IGelfand et al.l ( 2009 1. To minimize the number 
of free parameters in our model, we do not consider emis- 






























































In Section 13.11 we presen t the observed properties of 
G54.1+0.3, and in Section 13.21 describe the algorithm 
used to determining which combinations of model param¬ 
eters are able to reproduce them. Finally, we compare 
our results with similar work (Section 13.2.11) . 


3.1. Observed Properties 

G54.1+0.3 is one of the best stud ied PWNe in the 
Milky Way . Associated w ith radio (jGamilo et al.l (20021 ) 
and X-ray (|Lu et a.l .1120071) pulsar PSR J1930+1852, it is 
also detected across the electromagnetic spectrum. This 
PWN has a similar extent at both radio and X-ray en¬ 
ergies (|Lu et al.l 120011 : ILang et ahl I2010D . with a semi- 
major axis of ~ 1.'25 and a semi-minor axis of ~ l'O 
(jLaner et al.l 120101 ) . Since our model assumes a spheri¬ 
cally symmetric PWN (Section [2]), we set the angular 
size of the PWN 6 pwn our model must reproduce to the 
“average” of its measured semi-minor and semi-major 
axes, and use these to determine the 3er lower and upper 
limits on 9 pwn (Table [lj). We also re quire our model to 
reprod uce its volume-integ rated radio (ILang et ~aLll2010f) , 
X-rav (iTemim et al.ll2010n . and TeV O'-rav dAcciari et all 
l2010h properties, listed in Table [lj We do not at¬ 
tempt to reproduc e the mid-infrar e d (mid-IR) proper - 
ties of G54.1+0.3 (IKoo et, al.ll2f)M ITemim et al.l [201?)h 
since this emission is dominated by material shocked and 
heated by the expanding PWN. Because we are using 
a one-zone model (Section [2]), we also do not attempt 
to any reproduce spatial variation s in its emission (e.g., 
ILu et al.1 1200 lb ITemim et al.ll2010h . 

Lastly, we require our model to reproduce the size of 
the SNR. The SNR around PWN G54.1+0.3 has been 
detected at both rad io (|Lang et al.l 120101) and X-ray 
(IBocchino et al.l 120101) energies, each reporting a some¬ 
what different angular radius 9 snl . To resolve this dis¬ 
crepancy, we analyzed an archival D-array 1.4 GHz VLA 
observation of this PWN, estimating an SNR angular ra¬ 
dius of ~ 6'6. We then estimated the error on 9 snl by 
setting 3 a upper and lower limits to those r eported by 
ILang et all ([2010) and IBocchino et all (I2010f) . 

As listed in Table [T] our model has to reproduce twelve 
different observed quantities - equal to the number of 
model parameters. As a result, our fit has zero degrees 
of freedom. While the distance d to G54.1+0.3 is a 
free parameter in our model (Tabic [T]), the fitting al¬ 
gorithm described in Section I5~21 favors d = 4.5 — 9 kpc, 
as derived from an analysis of its Hi absorption spectrum 
(iLeahv et al.ll2008h . 


3.2. Model Fit 


( 4 ) 


l200l . we set the true age t age of G54.1+0.3 to: 

± _ 2T c h 

^age — T 7sd •> 

p -1 

and the initial spin-down luminosity Eq of this pulsar 
to: 

^0=^(1 + ^)^, ( 5 ) 

where p and r s d are respectively the pulsar’s braking in¬ 
dex and spin-down timescale. 

For a given combination, we first determine the model- 
predicted value of each observable M. We then calculate 
the likelihood C(V\9) this set of parameters accurately 
represents the data: 


am=U^r- 

V27T(J i 


1 ( Mi-TDi \ 5 

s 2 V ) 


( 6 ) 


where cq is the error on each observed quantity. 
The MCMC algorithm then searches the possible 12- 
dimensional parameter space for the combinations with 
the largest In C: 




In 


i=l 

1 

= -2 X " 


1 


■ 2 ' C, 

where y 2 is defined as: 
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i= 12 
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( 7 ) 

( 8 ) 

(9) 


and C is the same for all combinations. Therefore, max¬ 
imizing In C is equivalent to minimizing y 2 . It conducts 
this search using the following procedure: 


1. For a given combination 9 n , evaluate ln£„. 

2. Propose 9 n+1 which is 9 n + f(9), where f(9) is a set 
of random, zero-mean, Gaussian distributed num¬ 
bers whose width varies for each model parameter. 

3. Calculate ln£„ + i for the proposed 9 n+ 1 . 

£ 

4. If !f +1 > S, where S is a random, Gaussian dis- 

*-'71 

tributed number between 0 and 1, then 9 n +2 is cal¬ 
culated with 9 n+ 1 as a starting point. Otherwise, 
0„+2 is calculated with 9 n as a starting point. 

The width of f{9) was chosen such that the 9 n+ 1 satisfies 























































Table 1 

The Set of Model Parameters with the Lowest (x 2 ~ 4.10), Their Predicted Properties of G54.1+0.3, and the Observed Values of 

These Quantities. 


Model Input Parameters 


Predicted Observables 


Parameter 

Value 

Observable 

Observed Value 

Predicted Value 

log(E sn / 10 51 ergs) 

-0.03 

^snr 

616 ± 0.'4 

6.5 

log (Mej/M 0 ) 

1.34 

Ppwn 

1.'14±0.'04 

1.12 

l°g( n ism/Cm -3 ) 

-2.29 

1.4 GHz Flux Density 

433 ± 30 Jy 

429 Jy 

Distance (kpc) 

4.90 

4.7 GHz Flux Density 

327 ± 25 Jy 

329 Jy 

V 

2.94 

8.5 GHz Flux Density 

252 ± 20 Jy 

257 Jy 

logOWl year) 

2.90 

±X,2—10 

(5.43 ± 0.035) x IO " 12 -2^ 

5.43 x 10 ” 12 

logins) 

-3.14 

r 

2.09 ±0.01 

2.09 ±0.002 


log(£ min /GeV) 

l°g(£break/GeV) 

log(£ min /GeV) 

Pi 

P2 


1.05 311 GeV Photon Density 

3.45 492 GeV Photon Density 

6.98 780 GeV Photon Density 

1.84 1.2 TeV Photon Density 

2.77 3 TeV Photon Density 


(1.10 ±0.56) x 10~ n P i° to " 

v ' cm^ s TeV 

(4.2 ±1.4) X 10- 12 

v ' cm z .s TeV 

(1.12 ±0.45) x 10“ 12 

v ' cm z s TeV 

(6.2 ±1.7) x 10- 13 P S ot °i: P 

v ' cm; s TeV 

(3.9 ±2.1) x 10 -14 p %° to ™ 

v J cm^ s TeV 


0.80 x 10- 11 p +°° 3 

cm z s TeV 

3.1 x io- 12 p T to ™ 

cm* s TeV 

1.21 x 10- 12 p i oto “, 

cn+ s TeV 

4.9 x IO ” 13 ^ oto ” 3 w 

cm; s TeV 

7.2 x 10 ~ 14 p ^ ot °” 3 f 

cm^ s TeV 


Photon Energy [eV] 



Figure 1. The broadband spectral energy diagram of PWN 
G54.1+0.3 predicted by the model described in Section [2] for the 
parameters listed in Table [l] The red, green, and purple points 
are, respecti vely, the observed radio, X-ray, and TeV 7 -ray emis¬ 
sion (Section 13.1 1 Table [I])- 


in Figure [Q this set of parameters accurately reproduces 
the broadband SED of this PWN. 

Our search of parameter space allows us to estimate 
the (statistical) confidence interval of a given parame¬ 
ter by first ordering, from lowest to highest, its value 



Figure 2. The ISM density rii sm (color scale) for different values 
of the initial kinetic energy E sn and mass M e j of the supernova 
ejecta for trials with % 2 < 7.10 (the 3cr parameter space). 

tion coefficient r xy , defined to be: 

N 

E + - x){yi - y) 

i=1 , (io) 

J2(xi-x) 2 Jj2(yi-y) 2 

i=i y i=1 

between each pair of model parameters x and y , where x 
and y are their average values, Xi and yi are their values 
for a particular trial, and N is the number of trials, using 
only trials with y 2 < 7.10 (spanning the 3er parameter 




























Linear Pearson Correlation coefficient r X y 


Table 2 

(Equation {ToJ) between the model parameters as calculated for all trials with x 2 < 7.10. 



log E sn 

log Mej 

log ni sm 

P 

log r sd 

log VB 

log Emax 

log F min 

Pi 

log Ebreak 

P2 

d 

log E sn 

1.00 

0.87 

0.82 

- 0.12 

0.03 

0.22 

- 0.35 

0.13 

0.16 

0.18 

- 0.17 

0.23 

log Mej 

0.87 

1.00 

0.59 

- 0.18 

- 0.12 

- 0.07 

- 0.29 

- 0.19 

- 0.23 

- 0.15 

- 0.49 

0.37 

log n ism 

0.82 

0.59 

1.00 

- 0.31 

0.44 

- 0.01 

- 0.43 

0.13 

0.18 

0.04 

0.25 

- 0.33 

P 

- 0.12 

- 0.18 

- 0.31 

1.00 

-0.83 

- 0.04 

0.01 

0.01 

- 0.13 

- 0.01 

- 0.01 

0.09 

log r ad 

0.03 

- 0.12 

0.44 

-0.83 

1.00 

0.05 

- 0.11 

0.18 

0.36 

0.12 

0.42 

- 0.48 

log VB 

0.22 

- 0.07 

- 0.01 

- 0.04 

0.05 

1.00 

0.09 

0.49 

0.68 

0.73 

- 0.05 

0.54 

log E max 

- 0.35 

- 0.29 

- 0.43 

0.01 

- 0.11 

0.09 

1.00 

0.14 

0.04 

0.15 

0.13 

0.17 

log F m in 

0.13 

- 0.19 

0.13 

0.01 

0.18 

0.49 

0.14 

1.00 

0.78 

0.68 

0.35 

0.07 

Pi 

0.16 

- 0.23 

0.18 

- 0.13 

0.36 

0.68 

0.04 

0.78 

1.00 

0.93 

0.40 

0.12 

log E]-, reak 

0.18 

- 0.15 

0.04 

- 0.01 

0.12 

0.73 

0.15 

0.68 

0.93 

1.00 

0.28 

0.36 

P2 

- 0.17 

- 0.49 

0.25 

- 0.01 

0.42 

- 0.05 

0.13 

0.35 

0.40 

0.28 

1.00 

-0.71 

d 

0.23 

0.37 

- 0.33 

0.09 

- 0.48 

0.54 

0.17 

0.07 

0.12 

0.36 

-0.71 

1.00 


Note : Values in bold indicates that | r xy > 0.5|, indicating a significant degeneracy between the two parameters. 



Figure 3. The x 2 for trials with different values of the braking 
index p and spin-down timescale r s a of PSR J1930+1852, with red 
signifying a lower x 2 (better fit) and black a higher x 2 (worse fit). 
The dots indicate trials with x 2 < 5.10 and are included to better 
demonstrate the degeneracy between these two parameters. The 
dumpiness of these points primarily reflects the sampling of the 
parameter space by our MCMC algorithm. 

degenerate, with a more energetic supernova explosion 
requiring a larger ejecta mass occurring in a denser envi¬ 
ronment (Figure [2j). A similar degeneracy was reported 
in a recent analysis of Kes 75, wh ich discusses possibl e 
physical origins for this behavior (|Gelfand et al.|[20l4l) . 
The pulsar braking index p and spin-down timescale r s d 
are also strongly degenerate, with higher values of p re¬ 
quiring lower values of r s d (Figure [3j). The break energy 
Ehre&k in the spectrum of particles injected at the ter¬ 
mination shock strongly depends on the low energy par¬ 
ticle index pi, with higher values of F'break requiring a 
“softer” (higher values of pi) particle spectrum (Figure 


Figure 4. The x 2 f° r trials with different values of the break 
energy Ebreak and low-energy particle index pi , with red signifying 
a lower x 2 (better fit) and black a higher x 2 (worse fit). The 
dots indicate trials with x 2 < 5.10 and are included to better 
demonstrate the degeneracy between these two parameters. The 
dumpiness of these points primarily reflects the sampling of the 
parameter space by our MCMC algorithm. 
































Figure 6. The x 2 for trials with different values of the high- 
energy particle index p 2 and distance d to G54.1+0.3, with red 
signifying a lower x 2 (better fit) and black a higher x 2 (worse fit). 
The dots indicate trials with x 2 < 5.10 and are included to better 
demonstrate the degeneracy between these two parameters. The 
dumpiness of these points primarily reflects the sampling of the 
parameter space by our MCMC algorithm. 

energy -Ebreak (and correspondingly higher values of p \) 
and a larger distance d (and correspondingly lower values 
of P2 ) • 


3.2.1. Comparison with Other models 

In this section, we compare our results with those ob¬ 
tained using other models for the evolution of a PWN 
inside a SNR to determine how our analysis is affected by 
the assumptions made by our model described in Section 
[2] allowing us to estimate the systematic uncertainty of 
this approach. The results of these different models are 
p rovided in Table [3] 

iChevalieil (1200 Si) uses the measured spectral proper¬ 
ties and radius of this PWN and the spin-down proper¬ 
ties of the central pulsar to primarily estimate the birth 
properties of the neutron star, assuming M e j = 5 M 0 
and E sn = 10 51 ergs - a combination not favored by 
our fits (Figure [2]). He did not attempt to reproduce 
the broadband SED, and set p\ and P 2 to values in¬ 
ferred from single power-law fits to the observed radio 
and X-ray spectrum. While the value of P 2 derived from 
this method agrees with our value, the value of p± does 
not since, for E mm m 10 GeV, the SED predicted by 
our model contains a spectral break between 1.4 and 4.8 
GHz (Figure [I]). Additionally, he assumes that t]b = f 
(lChevalieill2005r ). significantly higher value than allowed 
by our fits. The higher value of t/b decreases the particle 
energy inside the PWN, resulting in an initial period Pq 


(0.5 — 1.6) x 10 51 ergs and t age ~ 1800 — 2400 years if 
the ions are ~ 2x hotter than the electrons. Both sets 
of E sn and f age are consistent with our results (Table 
[3]). They also found that M e j = 8 M 0 , p = 3, and 
r s d = 500 years can reproduce the radius of the PWN 
and SNR dBocchino et al.ll2010f) in agreement with our 
results. Since they did not attempt to reproduce the 
broadband SED of this source, this analysis does not 
constrain the magnetization or spectrum of particles in¬ 
jected into the PWN at the termina tion s h ock. _ 

G54 .1+0.3 was also analyzed bv iTanaka fe Takaharal 
(12011H . who reproduce both the size and broadband SED 
of this PWN using a model very similar to ours (Sec¬ 
tion E]) but include inverse Compton scattering of elec¬ 
trons off photon fields other than the CMB: an optical 
(T = 4000 K) photon field with an energy density it op t = 
0.5 eV cm' 3 , and an IR (T = 40 K) photon field with an 
energy density U\ r = 0.5 eV cm -3 or u; r = 2.0 eV cm' 3 

- finding that ? 7 b, Inbreak, and the parameters associated 
with the energetics of the n eutro n star (r sr i, t age , Eq, and 
Po) depend on u; r (jTanaka fc Takaharal 1201 il l. As listed 
in Table [U in general our parameters agree - though their 
analysis favors a higher value of Pq (less energetic neu¬ 
tron star) due to the inclusion of these additional photon 
fields. 

Similar results were obtained bv iTorres et all (2014), 
which uses an evolutionary model that includes the dif- 
fusion of particles both i nside and out of th e PWN 
(|Martin et al.l I2012H . Like ITanaka fe Takaharal (120111) , 
they include emission from electrons inverse Compton 
scattering off two photon fields in addition to the CMB, 
one with Xfir = 20 K and energy density rtpiR, = 
2.0 eV cm -3 and the other with T ntr = 3000 K and 
energy density ufir = 1-1 eV cm” 3 (ITorres et al.ll2014h 

- again deriving a lower Eq (higher Po) than our analysis. 
This model also assumes the maximum energy of parti¬ 
cles is limited by confinement in the termination shock - 
finding that the current value of P m ax is similar to what 
we require for our model. 

Lastly, we compare our results with those of iLi et ahl 
(120101 1. who model the broadband SED of G54.1+0.3 for 
both a leptonic and combined leptonic_and h adro nic ori¬ 
gin for the observed y-rays. Like ITorres et al.l (I2014H . 
they allowed the maximum energy of particles injected 
at the termination shock E max to vary, setting it to the 
e nergy w hose Larmor radius is the radius of the PWN 
(iLi et all 12010 1. Their model also allows leptons to es¬ 
cape from the PWN, and that these particles inverse 
Compton scatter off the CMB, background IR and op¬ 
tical photons from the Milky Way, and emission from 
the IR “lo op” and its embedded point sources around 
this PWN (iKoo et al.ll2008t iTemim et al.li20l6 ). In the 






















































Figure 7. The magnetization of the pulsar wind t]b (color scale) for different values of the break energy F^reak and low-energy particle 
index pi [right) and distance d and high energy particle index p 2 {left). Both are calculated for for trials with % 2 < 7.10 (the 3cr parameter 
space). 


Table 3 

The 90% Confidence Interval of the Properties of G54.1+0.3 Derived from our Analysis, Compared with Values Derived from Previous 

Analyses of this Source. 


Parameter 

This Work 

Chevalier (2005) 

Bocchino et al. (2010) 

Li et al. (2010) 

Tanaka Takahara (2011) 

Torres et al. (2014) 

Esn (10 51 ergs) 

0.08 - 1.5 

= 1 

0.3 - 1.6 



= 1 

"ej ("©) 

5.7 - 44 

= 5 

= 8 



= 20 

"ism C cm_3 ) 

(0.03 - 6.3) X 10 -3 


~ 0.2 



= 10 

Distance (kpc) 

4.6 - 8.1 

~ 5 

= 6.2 

= 6.2 

= 6.2 

= 6 

Braking Index p 

1.90 - 2.93 

= 3 

= 3 

= 3 

= 3 

= 3 

T sd (y ears ) 

280 - 3500 

« 1400 

= 3 


600 / 1200 

1171 


(0.44 - 2.2) X 10 -3 

= 1 


~ 1.5 X 10 -3 

0.3 X 10 -3 / 2 X 10 -3 

5 X 10 -3 

Emin (G«V) 

0.31 - 15 



= 0.05 

< 10 


Efareak ( T ®V) 

0.71 - 11 



= 0.26 

0.15 / 0.09 

0.3 

fimax (PeV) 

0.96 - 2700 



Variable 

> 0.5 

0.38 (Variable) 

PI 

1.43 - 2.08 

= 1.26 


= 1.2 

1.2 

1.2 

P2 

2.60 - 2.78 

= 2.8 


~ 2.8 

2.55 

2.8 

Age [years] 

2100 - 3600 

rs 1500 

1800 - 3300 

~ 2900 

2300 / 1700 

1700 

Eq [ergs s — - 1 ] 

(0.06 - 2.5) X 10 39 

« 5.1 X 10 37 

= 4 X 10 38 

= 1.4 X 10 39 

2.9 X 10 38 / 6.9 X 10 37 

7.2 X 10 37 

P 0 [ms] 

32 - 84 

w 100 

= 56 


62 / 87 

87 


Note : I Chevalier! d20051) do not specify a braking index p for this ne utron star, and the quoted values of r s d and Eq are calculated assuming 
p = 3 for the age derived in their analysis. As described in Section 13. 2. lllTanaka Sz Takaharal (120111) calculate the properties of this PWN 
assuming two different energy densities of the background IR photon field, with the values to the left of the “/” inferred for a lower energy 
density while the values to the right are those inferred for a higher energy density. 


the pulsar wind (Section l4~3l) . 

4.1. Progenitor Star 

The initial kinetic energy E sn and mass M e j ejected 
in a core-collapse supernova depends on the initial 
mass, metallicity, and evolution of the progenitor star 
(e.g., iHeeer et al.l 120031 ). G54.1+0.3 has a galactocen- 
tric radius (~ 6.5 — 7.5 kpc for the favored distance of 
d ~ 5 — 8 kpc) similar to the Sun’s (~ 8 — 8.5 kpc; 
lAndrievskv et al.ll2002alTbh . suggesting its progenitor had 
approximately Solar metallicity. A massive star in a bi- 


2. have extremely high (~ 50 Mq) initial mass but 
explodes after shedding much of this mass as a 
Wolf-Rayet star, resulting in a low (< 3 Mq) ejecta 
mass. 

As shown in Figure [2] a “canonical” supernova explosion 
energy of E sn ~ (0.3 — 1) x 10 51 ergs requires a higher 
ejecta mass (M e j > 10 M 0 ; Figure^). 

We can further constrain these parameters using the 
properties inferred from an analysis of the I R spe ctrum 
of the material su rrounding the PWN (|Koo et al.l 120081 

































Figure 8. The expected expansion speed of the ejecta just outside 
the PWN u e j(i? P wn) for different values of the initial kinetic energy 
E S n and mass M e j of the supernova ejecta for trials with % 2 < 7.10, 
M e j < 20 Mq, and p < 3. 

120031 ). This progenitor mass is further supported by the 
identification of O and B stars emb edded inside the SN 
ejecta dust surrounding this PWN (iTenfim et al.ll2010l) . 
Therefore, G54.1+0.3 was likely produced by the core¬ 
collapse of a ~ 15 — 20 Mq star in a massive star cluster 
possibly the most massive member of this cluster, and 
therefore the first to explode. 

This progenitor mass, and association with an mas¬ 
sive star cluster, can explain the low ISM density nj sm 
required by our model (Table 0. The winds of main- 
sequence massive stars are thought to creat e low-density 
bubble with a radius Pb (IChen et alJl2013l) : 

f -Pjsm/fcB \ 
\ 10 5 cm -3 K ) 


where P lsm is the pressure of the medium outside the 
wind bubble and fcs is Bo ltzm ann ’s con stant. For a dis¬ 
tance of ~ 4.5—9 kpc dLeahv et al .1120081) . this bubble will 
have an angular size of 0b ~ 217 — 14'. Winds from the 
additional massive stars in the cluster will only increase 
the size of this bubble, increasing the likelihood that the 
SNR is expanding inside a low density environment. 

4.2. Neutron Star Formation and Evolution 

The birth properties of a neutron star reflect the 
physics of its formation. The initial spin period Pq 
and surface magnetic field of the neutron star depend 


i?b = 


M 

(1.22 ± 0.04)— - (9.16 ± 1.77) 

Mq V ' 


'7—18 


-Pjsm/fcB 
10 5 cm -3 K 


pc, 


oo 


the neutron star is limited by gravitational waves re¬ 
sulting from r-rnode instabilities generated by the “fall¬ 
back” of material during the supernova onto the proto¬ 
neutron star, the B ns = 1.0 x 10 13 G dipole surface 
magnetic field strength inf erred from the_timing prop¬ 
erties of PSR .11930 +1852 (iGamilo et al.ll200l requires 
Po ~ 30 — 80 ms dWatts fe Anderssonll2002l) ' - consistent 
with the range favored by our model (Table [3]). 


4.3. Pulsar Wind 


The rotation of the neutron star generates a strong 
electric potential (voltage) <I> at its magnetic poles re¬ 
sponsible fo r bot h c reating part icles in its magnetosphere 
(e.g., IGoldreich fe Julian! I1969D . The pulsar wind con¬ 
sists of particles which exit the magnetosphere along 
open field lines, expected to occur at a minimum rate 

N gj 


1V G j = — = 7.6 x 10 
e 


33 


^45 


P 


P | 3 4 x 10 -13 s/s 


1;3) 


(IGoldreich fe Juhanlll969i : rBucciantini et al.l[201lD where 
the neutron star’s moment of inertia is I = I 45 x 
10 45 g cm 2 , P 33 = P/33 ms, and P is the neutron star’s 
period-derivative. However, how particles are both cre¬ 
ated and leave the neutron star magnetosphere is poorly 
understood. 

If particles are neither created nor destroyed between 
the light cylinder and the termination shock, we can cal¬ 
culate the rate particles leave the magnetosphere N for a 
particular trial using Equations [2] and [3] Our assumption 
that the parameters regulating the spectrum of particles 
injected at the termination shock (P m i n , Pbreak, Pmax, 
1 Pi, and P 2 ; Table [IJ are constant results in N ex E over 
3 . the life time of the PWN. As result, in our model the 
^ ^multiplicity of the pulsar wind k 


N 

+> ^ 5 + (14) 

varies with time. Therefore, in addition to calculating 
the current multiplicity Know, we also calcu l ate the time- 
integrated multiplicity Ki nt (e.g.- lde Jagerll2007l) : 


^age 

/ Ndt 

Kint = 7 -^-• (15) 

tage 

/ Nqj dt 
0 

Our analysis of G54.1+0.3 indicates that Know ~ 10 3 
«int ~ (1 — 3) x 10 5 - both in good agreement with 
the v alues deri ved f rom similar analyses of other PWNe 
































Near the neutron star, the pulsar wind is expected to 
be highly magnetized (jib ~ 1). However, our model 
requires that 773 ~ 10 -3 (Table [3J when the pulsar 
wind is injected into the PWN - requiring that mag¬ 
netic energy is converted to particle energy between 
the neutron star’s light cylinder and the termination 
shock (e.g., iKirk fc Skiasraasenll2003ll . Currently, mag¬ 
netic reconnection in this region is thought to trans¬ 
form the pulsar wind from a strongl y mag n etized to 
a weakly magnetized outflow (e . g.. IKirk fc Skiasraasenl 
1200. 'll : ISironi fc Snitkovskvl I2011L 1201 -111. Efficient mag- 
netic reconnection requires that (IKirk fc S kiacraascn 

IMl): 

/ 7T^6^ • \ 4 

^ < 3 ( — 2 - 5 E ) (16) 

\ TOg c° j 

where e and m e are, respectively, the charge and mass 
of a positron, c is the speed of light, and p,, the 
energy per un it mass en ergy of the pulsar wind, is 
(iKirk fc Skiaeraasenll2003lf : 


E 

Nmc 2 ’ 


(17) 


equivalent to the bulk Lorentz factor of the pulsar wind 
7 W before it reaches the termination shock (i.e., “up¬ 
stream” from the shock). Therefore, magnetic reconnec¬ 
tion is viable as long as the spin-down luminosity of PSR 
J1930+1852 is: 


E> 


ra); c 5 /u 4 
8l7r 3 e 2 


(0.2 - 1.3) x 10 35 —, (18) 

s 


for ji sa (1.5 — 2.5) x 10 5 as favored by our model. Since 
this critical E is well below its current E « 1.7 x 10 37 
(ICamilo et al.l2002[) , magnetic reconnection should occur 
in the pulsar wind before it reaches the termination shock 
- possibly explaining the weakly magnetized pulsar wind 
required by our model. 

Recent numerical simulations suggest that magnetic 
reconnection is the pulsar wind will produce particles 
whose spectrum is well described by a power-law with 
particle index p < 2, as required by our model for E < 
Kbreak (Table [3]), up to an energy : 


E, 


max,recon 


>m e c 


( Brecon + 1)(2 -p) 


ijP~ 1) 


(19) 


if the ratio of magnetic to particle energy in the magnetic 
recon nection region is a recon > 10 (ISironi fe Spitkovskvl 
I2014H . We can test if this is plausible calculating cr reC on 
if, in liquation 1191 p — Pi and E max recon Kbreak ■ 


iBucciantini et all[2011 ). where $ is the voltage of the 
pulsar’s magnetosphere: 


<f> = 



( 21 ) 


The cur rent sp in-down luminosity E of PSR J1930+1852 
(ICamilo et al Jl2002lt would suggest that E max ss 6 PeV 
in its magnetosphere - consistent with the values £ max 
required by our modeling (Table [3]). Another possibility 
is that these particles are created by additional accelera¬ 
tion at the termination shock. Simulations suggest that 
efficient acceleration of an electron- positro n plasma in 
this region requires r]B < 10~ 3 (e.g.. ISironi et al.ll2013h . 
again consistent with the range of values favored by our 
modeling. The maximum particle energy is expected be 
limited by either synchrotron cooling or diffusion away 
from the termination shock, with the theoretic maximum 
energy E m ax,theory being the lower of the two. For the 
pulsar wind properties favored by our modeling, the max¬ 
imum energy of the particles accelerated at the term ina¬ 
tion s hock is limited by diffusion, such that: (ISironi et al.l 

IM 1 : 


K maX: theory — 1-9 X 10 Ul e C 


E 


N 


!Q38.5 


1Q40 g—1 


( 22 ) 


H5-25 PeV. 


(23) 


Since ~ 50% of our trials have £ max < iJ maXi theory; our 
results are also consistent with highest energy particles 
being produced at the termination shock. 

Lastly, numerical simulations suggest the spectral 
shape of particles injected into the PWN at the termi¬ 
nation shock depends strongly on the structure of th e 
unshocked pulsar wind ('e.g.. ISironi &; Spitkovskvll2011lh 
When it leaves the neutron star magnetosphere, the pul¬ 
sar wind is expected to be primarily equatorial and com¬ 
posed o£n3giqns^ of alt ernating magnetic field directions 
(e.g., iBogovalovI 1 1 9991) of width A. The shape of the 
resul tant particle sp ectrum is expected to depend on 
(ISironi fe Spitkovskvll201ll ): 


-~ 47 TK 

TLU 


Rlc 
Rts ’ 


(24) 


where tl and a are, respectively, the relativistic Larmor 
radius and magnetization of the unshocked pulsar wind, w 
k is the multiplicity (Equation [14]) , i?xs is the radius of 
the termination shock, i?xs is th e radius of termination 
shock, and 7 ?lc is the radius of the light cylinder: 















































Ngj ~ 4.69 x 10 34 s -1 and Plc ~ 6.53 x 10 8 cm. Ad¬ 
ditionally, analysis of a Chandra observation identified a 
ring with semi-major axis $ts = 577 centered on the pul¬ 
sar, which is believed to m ark the position of the termi- 
nation shock in this PWN (iLu et al.ll2002t iTemim et al.l 
2010). For these values, the trial parameters with the 
lowest Y 2 favor ~ 10 -5 — 10 -4 , in contradiction with 

/v rr.cr 1 

the results of lSironi & Spitkovskvi ()2011[ ). 

5. OBSERVATIONAL TESTS 

While our evolutionary model for a PWN inside an 
SNR (Section [2]) reproduces the observed properties of 
G54.1+0.3 for a wide range in parameter space (Table 
H, it is important to test the validity of this model by 
predicting the value of additional observable properties. 
Thanks to our parameter exploration, not only can we 
predict the values of future observations, we can also es¬ 
timate the resulting improvement in the allowed physical 
parameters. For these predictions, we only use trials with 
X 2 < 7.10, M e j < 20 Mq, and p < 3.0. We only consider 
trials with M e j < 20M 0 since stellar evolution models 
suggests this is the ma ximum ejecta mass possible for a 
Solar metallicity star dWooslev et al.l l2002t iHeger et al.l 
2003; Heger 2015, private communication), and only tri¬ 
als with p < 3 since p > 3 has yet to be measure d from 
any isolated neutron star (e.g.. iLivingstonell20111 h 

Our model can predict properties of the SNR around 
G54.1+0.3 not yet measured, for example its expansion 
velocity u snr . Due to the young age and low ISM den¬ 
sity preferred by our model, we predict an extremely fast 
Tgnr ~ 3000 km s _1 - among th e highest measured or in - 
ferred for any other SNR (e.g.. [Ghavamian et al.l[2007th 
This suggests the identified radio and X-ray shell may 
not actually be a SNR but the progenitor’s stellar wind 
bubble (Section 14.11) . This can be determined by the 
measuring its radio spectral index (a, where flux density 
Sj, oc u a ), since the free-free emission expected to dom¬ 
inate the radio emission from a stellar wind bubble has 
a > 0 while SNRs typically have a ~ — 0.7. If future 
studies indicate this is a stellar wind bubble, our model 
would still favor a ~ 15 — 20 M 0 progenitor based on 
the properties of the IR emission around the PWN (Sec¬ 
tion |4T|), but would offer much weaker constraints on the 
density of the surrounding ISM. 

We can also predict currently unmeasured properties 
of the PWN, and determine what can be gained from 
their measurement. For example, our model predicts 
the average angular radius of the PWN is expanding 
by 0 P wn ~ 20 — 50 mas year -1 , and this value is sen¬ 
sitive to the mass M e j and initial kinetic energy E sn of 
the supernova ejecta because the PWN has not yet col¬ 
lided with the SNR reverse shock (Figure [9]). This is po- 



Figure 9. The expected angular expansion rate of this PWN’s 
radius 0 pw n for different values of the initial kinetic energy E s n 
and mass M e j of the supernova ejecta for trials with % 2 < 7.10, 
M e j < 20 Mq, and p < 3. 

tion shock E m \ n (Figure fill). All four of these quantities 
are m easurable by new observ ing facilities such as LO- 
FAR (Ivan Haarlem et iiD 201 Lastly, we find that the 
absorbed 5 — 80 keV f lux of G54.1+0.3, m easurable by 
the NuSTAR satellite (jHarrison et all 1201 31. is strongly 
depends with the distance to this source (Figure fl2l) 
likely a resu lt of the parameter degeneracies discussed in 
Section [3721 


6 . SUMMARY AND CONCLUSIONS 

In summary, we have fit the observed properties of 
G54.1+0.3 using a one-zone model for the evolution of a 
PWN inside an SNR (Section [2]). This model can repro¬ 
duce its observed properties (Section [3.1 1) , and suggests 
that the progenitor was an isolated ~ 15 M 0 star, most 
likely the member of a massive star cluster, which ex¬ 
ploded in a low density environment possibly produced 
by its stellar wind (Section 14.IF The resultant neu¬ 
tron star, PSR J1930+1852, had an initial spin period 
Po ~ 30 — 80 ms (Section PO!) . Our model requires that 
the current multiplicity of particle production in its mag¬ 
netosphere is k ~ (1 — 3) x 10 5 , and suggests that the 
magnetosphere electric potential is sufficient to acceler¬ 
ate particles to the highest energies E max required by 
our model. The low magnetization of the pulsar wind 
and low-energy component of particle spectrum can be 
attributed to acceleration resulting from magnetic recon¬ 
nection between the light cylinder and the termination 
shock, though our model suggests the “stripes” in the 
unshocked pulsar wind are too narrow for acceleration 
at the termination shock to nrnduce the broken newer- 
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Figure 10. The distance d to G54.1+0.3 (left) and initial spin period Pq of PSR J1930+1852 (right) for the values of the 60 MHz Sqq 
and 150 MHz S 150 flux densities predicted by trials with x 2 < 7.10, M e j < 20 Mq, and p < 3. 



Figure 11. The minimum energy in the pulsar wind E m i n for 
the 30—80 MHz ( 0 : 30 — 80 ) and 120—240 MHz (a 120 — 240 ) spectral 
indices of G54.1+0.3 predicted by trials with x 2 < 7.10, M e j < 
20 Mq, and p < 3. 

Aliu for information regarding the GeV spectrum, and 
Roger Chevalier, Vikram Dwarkadas, Daniel Patnaude, 
and Lorenzo Sironi for useful advice. 
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